/*
 *  lib-smv: Matrix - Vector Library targeted to Structural Mechanics Problems
 *  Copyright (C) 2006-2008  St.Brcic, Lj.Zugic
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
/** @file EigVal.cpp
 *  @brief Implementation of class EigVal - solution of the eigenvalue problem
 *
 *  EIGENVALUE PROBLEMS:  Standard and Generalized   
 *  A*x = L*x and A*x = L*B*x
 *
 *  Matrices A and B are real, square symmetric and positive definite (SPD).
 *
 *  Implemented routines are:
 *  (1) Jacobi transformation (only for standard problems)
 *  (2) Householder transformation combined with QL transformation
 *  (3) Lanczos algorithm (for standard and generalized problems) 
 */ 

#include "EigVal.h"

namespace smv_lib
{

//-------------------------------- constructors
//! empty (default) constructor
EigVal::EigVal(void)
{
  isMatrix_A = false;
  isMatrix_B = false;
  isNumber   = false;
  isStandard = false;
  isGeneral  = false;
  isMethod   = false;
  isNameGiven = false;
  isSolution = false;
}

//! constructor with one matrix and a number of eigenvalue/vector pairs
//! for STANDARD eigenvalue problem
//
EigVal::EigVal(Matrix& a, int n) :
    A(a), EV_Number(n)
{
  N = A.getRows();
  assert (N == A.getCols());  // Matrix must be square
  assert(EV_Number <= N);
  // all OK so far
  isMatrix_A = true;
  isNumber   = true;
  isMethod   = false;
  EV_Type = STANDARD;
  if (EV_Number == N)
    EV_Anal = COMPLETE;
  else
    EV_Anal = PARTIAL;
  isSolution = false;
}

//! constructor with two matrices and a number of eigenvalue/vector pairs
//! for GENERALIZED eigenvalue problem
//
EigVal::EigVal(Matrix& a, Matrix& b, int n) :
    A(a), B(b), EV_Number(n)
{
  N = A.getRows();
  assert (N == A.getCols());  // Matrix A must be square
  int m = B.getRows();
  assert (m == B.getCols());  // Matrix B must be square
  assert(N == m);             // Orders of A and B must be the same
  assert(EV_Number <= N);     // Number of eigenvalues
  // all OK so far
  isMatrix_A = true;
  isMatrix_B = true;
  isNumber   = true;
  isMethod   = false;
  EV_Type = GENERAL;
  if (EV_Number == N)
    EV_Anal = COMPLETE;
  else
    EV_Anal = PARTIAL;
  isSolution = false;
}

//-------------------------- destructor
//
EigVal::~EigVal(void)
{
  if (isSolution) {
    delete [] xi;
    xi = 0;
    delete m_Jacobi;
  	m_Jacobi = 0;
	delete    m_TQLI; 
    m_TQLI = 0; 
    delete m_Lanczos;
    m_Lanczos = 0;    
  }
}

//! set Matrix A (for standard problem)
//
void EigVal::setMatrix_A(Matrix& a)
{
  A = a;
  N = a.getRows();
  assert(N == a.getCols());
  isMatrix_A = true;
  EV_Type = STANDARD;
  isSolution = false;
}

//! set Matrices A and B (for generalized problem)
//
void EigVal::setMatrix_B(Matrix& a, Matrix& b)
{
  A = a;
  B = b;
  N = A.getRows();
  assert(N == A.getCols());
  int m = B.getRows();
  assert(m == B.getCols());
  assert(N == m);
  isMatrix_A = true;
  isMatrix_B = true;
  EV_Type = GENERAL;
  isSolution = false;
}

//! set number of eigenvalue/vector pairs
//
void EigVal::setNumber(int num)
{
  EV_Number = num;
  assert (EV_Number <= N);
  isNumber   = true;
  isSolution = false;
}

//! set the name of the eigenvalue problem
void EigVal::setName(string name)
{
  Name = name;
  isNameGiven = true;
}

//=================================== SOLUTION OF THE PROBLEM
//
//! solve the eigenvalue problem and write results to the file
void EigVal::solveEV_Problem(string fileName)
{
  // check input data
  if (checkInput() == false) {
     cerr << "\nEigVal: Problem with INPUT data" << endl;
     return;
  }
  // select method of solution
  selectMethod();
  // now, obtain the solution
  isSolution = true;
  // write results to the file
  writeEV_Results(fileName);
}

//! check input data
//
bool EigVal::checkInput()
{
  bool test = false;
  switch (EV_Type) {
      case STANDARD:
         if (isMatrix_A && isStandard) {
           test = true;
           isMethod = true;
         }
         break;
      case GENERAL:
         if (isMatrix_A && isMatrix_B && isGeneral) {
           test = true;
           isMethod = true;
         }
         break;
  }
  if (isNumber == false || isMethod == false)
     return test;
  return test;
}

//! select the method of solution
//
void EigVal::selectMethod()
{
   // select the method of solution
   switch (EV_Type) {
      case STANDARD:
         switch (EV_Standard) {
            case sJACOBI:  standardJacobi();
                           break;
            case HOUSE:    Householder_TQLI();
                           break;
            case sLANCZOS: standardLanczos();
                           break;
         }
         break;
      case GENERAL:
         switch (EV_General) {
            case eqSTAND:  equiStandard();
                           break;
            case gJACOBI:  generalJacobi();
                           break;
            case gLANCZOS: generalLanczos();
                           break;
         }
         break;
   }
}

//--------------------------------- various methods for eigenvalue solution
//! standard Jacobi transformation
void EigVal::standardJacobi()
{
   m_Jacobi = new Jacobi(A);
   m_Jacobi->calculate();
   X   = m_Jacobi->giveModalMatrix();
   Lam = m_Jacobi->giveSpectrum();
}

//! Householder reduction with QL transformation
void EigVal::Householder_TQLI()
{
   m_TQLI = new TQLI(A);
   m_TQLI->calculate();
   X   = m_TQLI->giveModalMatrix();
   Lam = m_TQLI->giveSpectrum();
}

//! standard Lanczos method
void EigVal::standardLanczos()
{
   m_Lanczos = new Lanczos(A,EV_Number);
   m_Lanczos->calculate();
   X   = m_Lanczos->giveModalMatrix();
   Lam = m_Lanczos->giveSpectrum();
}

//! transformation to equivalent standard problem: (B_inv * A - Lambda I)x = 0
void EigVal::equiStandard()
{
   Cholesky* m_Chol= new Cholesky(B,chol_STANDARD); 
   assert (m_Chol!= NULL);
   m_Chol->Decomposition();
   m_Chol->Inverse();
   Matrix L_inv = m_Chol->get_L_inv();
   L_inv.set_MatType(GEN_TYPE);
       
   Matrix Temp = Dot(L_inv,A);
   Matrix Temp1 = Tra(L_inv);
   Matrix C =  Dot(Temp,Temp1);
   C.set_MatType(SYMMETRIC); // C is always symmetric
        
   m_Jacobi = new Jacobi(C);
   m_Jacobi->calculate();   
   Lam = m_Jacobi->giveSpectrum();
   Matrix LTX   = m_Jacobi->giveModalMatrix(); // Tra(L)*X
   
   X = Dot(Temp1,LTX );       
   delete m_Chol;
   m_Chol = 0;                
}

//! generalized Jacobi transformation
void EigVal::generalJacobi()
{
   m_Jacobi = new Jacobi(A,B);
   m_Jacobi->calculate();
   X   = m_Jacobi->giveModalMatrix();
   Lam = m_Jacobi->giveSpectrum();
}

//! generalized Lanczos method
void EigVal::generalLanczos()
{
   m_Lanczos = new Lanczos(A,B,EV_Number);
   m_Lanczos->calculate();
   X   = m_Lanczos->giveModalMatrix();
   Lam = m_Lanczos->giveSpectrum();
}

//================================= SOLUTION RESULTS
//! reluts are given in the file
// (private function)
void EigVal::writeEV_Results(string fileName)
{
  
  string FileName, message;
  if (fileName != "")           // if the output file name is explicitely given
     FileName = fileName;
  else {
  if (isNameGiven == true)      // otherwise,
     FileName = Name + ".txt";  // output file name is according to the problem name
  else
    FileName = "EigVal.txt";    // dafault output file name
  }
  //--- output file
  ofstream out;
  out.open(FileName.c_str(),ios::out);
  out << "================== EIGENVALUE ANALYSIS ========================= \n\n";
  if (isNameGiven)
     out << "           Problem name      = " << Name << endl << endl;
  else
     out << "           Problem name      = " << "Eigenvalue solution"  << endl << endl;
  out << "           Basic parameters: " << endl;
  switch (EV_Anal) {
    case COMPLETE:
      message = "Complete eigenvalue problem";
      break;
    case PARTIAL:
      message = "Partial eigenvalue problem";
      break;
  }
  out << "(1) --- " << message << endl;
  switch (EV_Type) {
    case STANDARD:
      message = "Standard eigenvalue problem";
      break;
    case GENERAL:
      message = "Generalized eigenvalue problem";
      break;
  }
  out << "(2) --- " << message << endl;
  if (EV_Type == STANDARD) {
    switch (EV_Standard) {
      case sJACOBI:
        message = "Jacobi transformation";
        break;
      case HOUSE:
        message = "Householder reduction with QL transformation";
        break;
      case sLANCZOS:
        message = "Lanczos metrhod";
        break;
    }
  } else if (EV_Type == GENERAL) {
    switch (EV_General) {
      case eqSTAND:
        message = "Transformation to equivalent standard problem";
        break;
      case gJACOBI:
        message = "Generalized Jacobi transformation";
        break;
      case gLANCZOS:
        message = "Generalized Lanczos method";
        break;
    }
  }
  out << "(3) --- " << message << endl << endl;
  out << "\n---- Results are: " << endl;
  out << "\n(1) Eigenvalues are: \n";
  out << Lam << endl;       // maybe better: out << getEV_Spectar() << endl;
  out << "\n(2) Eigenvector matrix is: \n";
  out << X << endl;         // maybe better: out << getEV_Matrix() << endl;
  out.close();
}

//! complete modal matrix is given
//
Matrix EigVal::getEV_Matrix()
{
  if (isSolution == true) {
    Matrix ans(X);
    return (ans);
  }
  else {
    Matrix ans(N, EV_Number); // empty matrix
    return (ans);
  }
}

//! complete eigenvalue spectar is given
//
Vector EigVal::getEV_Spectar()
{
  if (isSolution == true) {
    Vector ans(Lam);
    return (ans);
  }
  else {
    Vector ans(EV_Number);  // empty vector
    return (ans);
  }
}

//! only one eigenvector is given
//
Vector EigVal::getEV_Vector(int i)
{
  assert(i>= 0 && i < EV_Number);
  Vector ans(N);
  if (isSolution == true)
    ans = X.getColData(i);
  return (ans);
}

//! only one eigenvalue is given
//
double EigVal::getEV_Value(int i)
{
  assert(i>= 0 && i < EV_Number);
  double ans = 0.0;
  if (isSolution == true)
    ans = Lam.getData(i);
  return ans;
}

//======================================================== NOT YET
//
///*------------- elmhes                       */
//
//#define SWAP(g,h)  {y=(g);(g)=(h);(h)=y;}
//
//void EigVal::elmhes (double **a, long n)
//
///* Reduction to Hessenberg form by the elimination method.
//   The real, nonsymmetric matrix a[1..n][1..n] is replaced by
//   an upper Hessenberg matrix with identical eigenvalues.
//   On output, thr Hessenberg matrix is in elements a[i,j], i<=j+1.
//*/
//{
//   long m, j, i;
//   double y, x;
//
//   for (m=2; m<n; m++)  {
//       x = 0.0;
//       i = m;
//       for (j=m; j<=n; j++)  {
//         if (fabs(a[j][m-1]) > fabs(x))  {
//            x = a[j][m-1];
//            i = j;
//         }
//       }
//       if (i != m)  {
//        for (j=m-1; j<=n; j++)   SWAP(a[i][j],a[m][j])
//        for (j=1; j<=n; j++)     SWAP(a[j][i],a[j][m])
//       }
//       if (x)  {
//        for (i=m+1; i<=n; i++)  {
//            if (y = a[i][m-1])  {
//               y /= x;
//               a[i][m-1] = y;
//               for (j=m; j<=n; j++)
//                   a[i][j] -= y * a[m][j];
//               for (j=1; j<=n; j++)
//                   a[j][m] += y * a[j][i];
//            }
//        }
//       }
//   }
//}
//
////------------- hqr
//
//#define SIGN(a,b) ((b) > 0 ? fabs(a) : -fabs(a))
//
//void EigVal::hqr (double **a, long n, double wr[], double wi[])
//
///* QR algorithm of an upper Hessenberg matrix a[1..n][1..n]
//   - only eigenvalues are determined.
//   On input a can be exactly as output from "elmhes".
//   The real and imaginary parts of the eigenvalues are returned
//   in wr[1..n] anr wi[1..n] respectively.
//*/
//{
//   long nn, m, l, k, j, its, i, mmin;
//   double z, y, x, w, v, u, t, s, r, q, p, anorm;
//
//   anorm = fabs(a[1][1]);
//   for (i=2; i<=n; i++)
//       for (j=(i-1); j<=n; j++)
//         anorm += fabs(a[i][j]);
//   nn = n;
//   t = 0.0;
//   while (nn >= 1)  {
//      its = 0;
//      do  {
//        for (l=nn; l>=2; l--)  {
//            s = fabs(a[l-1][l-1]) + fabs(a[l][l]);
//            if (s == 0.0)  s = anorm;
//            if ((double)(fabs(a[l][l-1]) + s) == s)  break;
//        }
//        x = a[nn][nn];
//        if (l == nn)  {
//           wr[nn]   = x + t;
//           wi[nn--] = 0.0;
//        }  else  {
//           y = a[nn-1][nn-1];
//           w = a[nn][nn-1] * a[nn-1][nn];
//           if (l == (nn-1))  {
//              p = 0.5 * (y-x);
//              q = p * p + w;
//              z = sqrt(fabs(q));
//              x += t;
//              if (q >= 0.0)  {
//                 z = p + SIGN(z,p);
//                 wr[nn-1] = wr[nn] = x + z;
//                 if (z)  wr[nn] = x - w/z;
//                 wi[nn-1] = wi[nn] = 0.0;
//              }  else  {
//                 wr[nn-1] = wr[nn] = x + p;
//                 wi[nn-1] = -(wi[nn] = z);
//              }
//              nn -= 2;
//           }  else  {
//              if (its == 30) util.nrerror("Too many iterations in HQR");
//              if (its == 10 || its == 20)  {
//                 t += x;
//                 for (i=1; i<=nn; i++)  a[i][i] -= x;
//                 s = fabs(a[nn][nn-1]) + fabs(a[nn-1][nn-2]);
//                 y = x = 0.75 * s;
//                 w = -0.4375 * s * s;
//              }
//              ++its;
//              for (m=(nn-2); m>=l; m--)  {
//                  z = a[m][m];
//                  r = x - z;
//                  s = y - z;
//                  p = (r*s-w)/a[m+1][m] + a[m][m+1];
//                  q = a[m+1][m+1] - z - r - s;
//                  r = a[m+2][m+1];
//                  s = fabs(p) + fabs(q) + fabs(r);
//                  p /= s;
//                  q /= s;
//                  r /= s;
//                  if (m == l)  break;
//                  u = fabs(a[m][m-1]) * (fabs(q) + fabs(r));
//                  v = fabs(p)*(fabs(a[m-1][m-1])+fabs(z)+fabs(a[m+1][m+1]));
//                  if ((double)(u + v) == v)  break;
//              }
//              for (i=m+2; i<=nn; i++)  {
//                  a[i][i-2] = 0.0;
//                  if (i != (m+2))  a[i][i-3] = 0.0;
//              }
//              for (k=m; k<=nn-1; k++)  {
//                  if (k != m)  {
//                     p = a[k][k-1];
//                     q = a[k+1][k-1];
//                     r = 0.0;
//                     if (k != (nn-1))  r = a[k+2][k-1];
//                     if (x=fabs(p)+fabs(q)+fabs(r))  {
//                        p /= x;
//                        q /= x;
//                        r /= x;
//                     }
//                  }
//                  if (s=SIGN(sqrt(p*p+q*q+r*r),p))  {
//                     if (k == m)  {
//                        if (l != m)
//                           a[k][k-1] = -a[k][k-1];
//                     }  else
//                           a[k][k-1] = -s * x;
//                     p += s;
//                     x  = p / s;
//                     y  = q / s;
//                     z  = r / s;
//                     q /= p;
//                     r /= p;
//                     for (j=k; j<=nn; j++)  {
//                         p = a[k][j] + q*a[k+1][j];
//                         if (k != (nn-1))  {
//                            p += r * a[k+2][j];
//                            a[k+2][j] -= p * z;
//                         }
//                         a[k+1][j] -= p * y;
//                         a[k][j]   -= p * x;
//                     }
//                     mmin = nn<k+3 ? nn : k+3;
//                     for (i=l; i<=mmin; i++)  {
//                         p=x*a[i][k] + y*a[i][k+1];
//                         if (k != (nn-1))  {
//                            p += z * a[i][k+2];
//                            a[i][k+2] -= p * r;
//                         }
//                         a[i][k+1] -= p * q;
//                         a[i][k]   -= p;
//                     }
//                  }
//              }
//           }
//        }
//      }   while (l < nn-1);
//   }
//}
//

}   // end of namespace smv_lib

//============================== end of file
